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Abstract. In this paper we develop two mathematical models to predict the 
^ release kinetics of a water soluble drug from a polymer/excipient matrix tablet. 

The first of our models consists of a random walk on a weighted graph, where 
the vertices of the graph represent particles of drug, excipient and polymer, 
respectively. The graph itself is the contact graph of a multidispcrsc random 
sphere packing. The second model describes the dissolution and the subsequent 
diffusion of the active drug out of a porous matrix using a system of partial 
I differential equations. The predictions of both models show good qualitative 

qq agreement with experimental release curves. The models will provide tools for 

designing better controlled release devices. 

• • 

. ^ 1. Introduction. The first Workshop on the Application of Mathematics to Prob- 

lems in Biomedicine took place from December 17-19, 2007 at the University of 
^ Otago in Dunedin, New Zealand. During that workshop, our group of nine, com- 

prising pharmaceutical scientists and mathematicians worked on the problem of 
predicting the drug release kinetics from a sustained release matrix tablet. The 
present paper grew out of this event. 

Sustained release (also called extended release) tablets are a common dosage 
form. A sustained release (SR) tablet is typically designed to release drug over 
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12-24 hours and might contain three times the dose of drug that is contained in an 
immediate release tablet. In this way, a patient need take a tablet only once a day 
rather than three times a day if immediate release tablets were used. This not only 
has the advantage of convenience for the patient but ideally provides more constant 
levels of drug in the body. Fluctuating drug levels can result in the patient being 
exposed to levels of drug which are too high at times, leading to harmful side-effects, 
and sub-therapeutic levels at other times. Sustained release tablets can smooth 
these fluctuations leading to better control of the patient's illness or symptoms. 

Drugs formulated in sustained release tablets typically have half-lives in the body 
of 3-8 hours. A drug with a long half-life given as a standard tablet once or twice a 
day inherently maintains reasonably constant levels in the body making sustained 
release tablets unnecessary. A drug with a very short half-life (say 1 h) would have 
to be given many times each day if formulated as a standard tablet. This might 
then suggest that an sustained release tablet would be useful, but such a tablet 
would have to contain perhaps 12 standard doses, and this would be potentially 
dangerous. 

Sustained release tablets can be formulated in various ways but the preferred 
approach, because of its apparent simplicity and ease of manufacture, is to make a 
sustained release matrix tablet. A matrix tablet can be made by mixing the drug 
with suitable excipients (non-active components of the formulation) and compress- 
ing the mix in a die at high pressure, say 100 MPa, thereby producing a tablet. 
Excipients are added to tablets for various reasons. Some are diluents to increase 
bulk and aid compaction. Others help in manufacture, for example lubricants and 
flow enhancers, while others influence the behavior of the tablet in water. For ex- 
ample, disintegrants added to standard tablets making them break-up when placed 
in water, are not used in sustained release tablets designed to remain intact. Ex- 
cipients vary in their physical properties, such as water solubility and mechanical 
properties. 

An important excipient in sustained release matrix tablets is the agent to control 
the drug release. Broadly, either lipids, such as fats and waxes, or polymers are 
used. In this paper we will further divide excipients into the insoluble polymer 
which forms the coherent matrix and soluble excipients. Some polymers, such as 
hydrogels, swell on exposure to water and create a rate-controlling gelatinous layer 
around the SR tablet. Other polymers, as used in the tablets described below, 
do not swell markedly but create a matrix through which water can penetrate to 
dissolve the drug particles and any water-soluble excipient particles. The dissolved 
particles then exist as molecules which diffuse down the concentration gradient, 
through the water-filled pores, to be released from the tablet. Several of the authors 
have been studying matrix tablets made by mixing drug with a polymer and one 
additional excipient, either lactose, a water-soluble brittle excipient whose particles 
might fracture during compression, or mannitol, a water-soluble plastic excipient 
whose particles plastically deform during compression. Other factors can potentially 
affect certain characteristics such as the hardness or sustained release behavior of the 
tablet. These include the particle size of the excipients, proportions of ingredients, 
compression force and post-compression thermal treatment. The latter factor can 
change the table porosity due to coalescence of polymer particles to form a coherent 
polymer network [1,8, 15]. The amount of drug in the tablet is determined by the 
total dose to be given to the patient, however, the proportions of the excipients vary 
with their function and the specific: tablet formulation. For example, a lubricant 
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such as magnesium stearate might be used at a 1% level, whereas the level of a 
rate-controlling polymer would be much higher (10-50%) but preferably the level 
would be kept as low as possible for cost reasons. 

When such tablets are placed in water, or in the intestinal fluid of the patient, the 
coherent polymer matrix sustains release of the drug by keeping the tablet largely 
intact and by providing a tortuous network through which water penetrates and 
dissolved drug and soluble excipient molecules diffuse out. That is, the postulated 
release mechanism involves penetration of fluid, dissolution of the drug and soluble 
excipient in the fluid, and outward diffusion of molecules of dissolved drug due to 
the concentration difference between the solution in the tablet and the intestinal 
fluid. Once released from the tablet, the drug is rapidly absorbed through the 
patient's intestine into the blood stream. Ideally, the delivery rate of the drug from 
the tablet should be such as to maintain the preferred blood level for an extended 
time. The left panel of figure 1 shows the structure of the polymer matrix of a 
tablet after dissolution of the solvable excipient and the drug. 

In this paper we propose two mathematical models that predict the release ki- 
netics of matrix tablets, taking into account parameters such as the composition 
of the powder mixture, powder particle sizes and applied curing temperature. The 
first model models the diffusion process as a random walk through the channels in 
the tablet that are formed in the compaction and heating process. In this model, 
the particle distribution within the tablet is represented by a graph embedded in 
R 3 whose vertices are generated from a random sphere packing and in which two 
vertices are joined by an edge if the corresponding spheres are in close proximity. 
On this contact graph we perform random walks that start at the location of each 
drug particle. These random walks mimic the diffusion of the drug molecules from 
their initial position to the edge of the tablet. By adjusting the probabilities of 
moving to vertices representing different types of particles, the composition of the 
tablet affects the length distribution of these exit paths in a realistic way. The right 
panel of figure 1 shows a sample path through a random packing of disks. 




FIGURE 1 . Left panel: Field emission scanning electron microscope 
(SEM) image of a matrix tablet after dissolution of drug and excip- 
ient at 400 x magnification. Right panel: A sample path through 
a schematic two-dimensional tablet from an inner particle to the 
edge. Small solid circles represent polymer particles. Notice that 
there are multiple paths from each initial particle location. 
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Low levels of drug in some polymer systems can lead to trapping and incomplete 
release. Percolation theory has been applied to understand these processes [3-6,10, 
17]. Monte Carlo simulations of drug release from a binary powder mixtures with 
inert (water insoluble) excipient have been performed by Villalobos et al. [30,31]. 
The drug and excipient particles were placed on the sites of a cubic lattice and 
the drug particles performed random walks on the accessible empty sites under 
a volume exclusion constraint. The authors reported different exponents for the 
release profiles as function of time, with v^-like behavior for large drug fractions. 
In our work we take into account drug, excipient and polymer, where the polymer 
plays the role of the inert matrix. More importantly, we also allow for differences 
in particle sizes and irregular placements of the particles in space. 

The second model investigates a different approach to predicting release kinetics 
of tablets using partial differential equations, where time and space are treated as 
continuous variables. Previous work using this approach includes [16,24]. Suitable 
versions of the diffusion equation (often in cylindrical coordinates, resembling the 
natural shape of a tablet) are formulated for the diffusion of dissolved drug and 
(e.g. in [24]) water penetrating the tablet. These models can take into account the 
dissolution kinetics of the drug in different types of buffer. Siepmann and Peppas 
[24] consider hydrophilic polymer matrices that swell and then dissolve, unlike the 
polymers used in the tablets that we are studying. While these results may not be 
directly comparable to our results for inert matrices, they do show that an increasing 
drug load results in faster release kinetics [24, Fig. 4]. 



2. Model 1: The diffusion of drug molecules is a random walk on a 
weighted graph. This first model is constructed in three distinct steps. In the 
first step, the positions of the particles in the tablet are taken from a random sphere 
packing, where each particle is a sphere of a fixed radius. In Step 2, this random 
sphere packing is represented as a graph embedded in K 3 , whose vertices are the 
particles, with edges between particles that are close to each other in space. The 
compression of the tablet is modeled as a deformation of this graph, and the heating 
of the tablet is modeled by removing edges. Finally, in Step 3, the diffusion of the 
drug particles to the exterior of the tablet is modeled as a random walk on the graph 
generated in Step 2. Drug release profiles are then generated as the distribution of 
exit times of these particles, using the following argument. 

Since we model the diffusion process as a random walk, the distance that a 
particular drug molecule must traverse from its position inside the tablet to the 
exterior is a random variable, L. Let /l(s) denote the density of these path lengths. 
Assume that a dissolved drug molecule performs Brownian motion with diffusion 
constant D along these paths. If the particle has to travel along a path of length 
s to get to the exterior of the tablet, then the travel time T is a random variable 
with density function (see [11, VI. 2]) 

S -a 2 

h\s{r) = /n ^_ 3 e^- 



V2ttDt 



In other words, Jt\ s is the density of T conditioned on the path-length being equal 
to s. Thus, we can integrate this conditional density against the density of path 
lengths to get the density of escape times of the particles from the tablet 

/T(r) = y o fT\ s (r)f L (s)ds= -—== e -^f L {a) da. 
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It's cumulative distribution function 

F T {t)= f / T (r)dr 
Jo 

gives the release profile of the matrix tablet. In Step 3, we estimate Fx(t) by 
performing many random walks on the graph created in Step 2. 

We now describe each of these three steps in more detail. Results from prelim- 
inary implementations of this model are given in Section 4. In this work we are 
not interested in supplying the optimal algorithm for completing each step in the 
modeling process. Rather, we are trying to validate the idea and, if the results 
appear sound, we can then optimize the algorithms. We assume that drug, excip- 
ient and polymer particles can be modeled as perfect hard (i.e. non-deformable) 
spheres, where all particles of each type have the same radius, but the three radii of 
the three different types of particles may differ from one another. The assumption 
that all particles of one species have the same radius is a simplification. In prac- 
tice, their sizes are distributed around a mean with a given variance. Relaxing this 
assumption, as well as including particles of other, non-spherical shapes are topics 
of future research. 

Step 1: We assume that prior to compaction the powder mixture is a dense 
multidisperse random sphere packing. Dense packings of disks or spheres have at- 
tracted the interest of mathematicians, physicists and engineers for a long time, see 
e.g. [9,14,18,27,29] and the references therein. It turns out to be extremely difficult 
to give a precise mathematical meaning to the concept of a "dense random sphere 
packing" , as the desired properties of high density and high degree of randomness 
conflict with each other [29]. The more dense a packing is, the more ordered it 
tends to be, culminating in the highly ordered lattice packings such as the tridiago- 
nal packing of disks in M 2 and the face-centered cubic packing of balls in K 3 . (The 
optimality of the latter packing was conjectured by Johannes Kepler in 1611 and 
proved by Thomas Hales in 1998.) While this question is certainly interesting from 
a mathematical point of view, we shall not pursue it further. Rather, we will work 
with the outcome of an experimental protocol that produces jammed configurations 
of hard spheres. 

The protocol used in this paper was first suggested by Lubachevsky and Stillingcr 
[18] and later expanded by Knott and coworkers [14]. The paper by Knott et 
al. [14] discusses the problem of packing spheres of oxidizer and fuel in a solid rocket 
propellant. Briefly, the Lubachevsky-Stillinger (LS) protocol proceeds as follows. 
A random initial configuration of sphere centers (xi(0), ^(O), ... xn(0)) and a 
random set of initial velocities (ui(0), i>2(0), . . . vjv(O)) are chosen. The spheres 
are partitioned into M radius classes, each class contains Afj spheres (i — 1, . . . , M) 
and N — JZj—i Nj. The spheres move within a container and collide with each 
other and with the walls of the container. As time evolves, the radius of class i 
spheres grows according to ri(t) = a»t. Thus the ratios ri(t)/rj(t) are preserved 
throughout the process. An event is a collision of two spheres or the collision of a 
sphere with the wall. At each collision of two spheres the velocities of the colliding 
spheres are updated in such a way that the spheres move away from each other 
after the collision. This requires an increase in mechanical energy since the radii 
are growing at a constant rate. Eventually, as the packing becomes more and more 
dense, the time between two collisions approaches zero and the process is stopped. 
The resulting configuration of spheres is (nearly) rigid, in the sense that no sphere 
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can be moved while keeping the centers of all other spheres fixed. For further 
details we refer to [14,18]. It should be noted that there exist alternative methods 
to construct random dense sphere packings [13,25,28]. 




Figure 2. Left panel: The result of the Lubachevsky-Stillinger 
protocol using a total of 735 spheres where drug, polymer and 
excipient particles are present in a mixture 90 : 100 : 545. The 
excipient particles are slightly larger than the others (see table 
2). The packing density is approximately 0.54. Right panel: The 
result of the Lubachevsky-Stillinger protocol using a total of 198 
spheres where drug, polymer and excipient particles are present 
in a mixture 90 : 100 : 8. The excipient particles are about 4 
times bigger than the others (see table 2). The packing density is 
approximately 0.63. 

Step 2: Having obtained a multidisperse dense random packing of spheres, this 
sphere packing is compressed as the powder is compressed in the die [25]. We 
model the die as a right circular cylinder whose axis is the z-axis and whose lower 
surface lies and remains in the xy-plane. The top surface of the cylinder which 
was originally in the plane {z — 1} is displaced by a < 1. A particle originally 
located at the level z is displaced to the level az. Using the new coordinates, a 
graph G is constructed, the proximity graph of the compacted sphere packing. This 
is a Euclidean graph G = {V,E} whose vertices are the compressed sphere centers. 
Each vertex carries a label L : V — > {D, P, X} that indicates whether the sphere is 
a drug, polymer or excipient particle, respectively. Two vertices are connected by 
an edge if the corresponding spheres are sufficiently close to each other. We also 
introduce a surface indicator s : V — > {0, 1} where s(v) = 1 if the sphere centered 
at v is in contact with the exterior of the tablet. 

Next, we form the heated contact graph G. If the centers of two spheres are 
connected by an edge, this edge is removed with a probability p that is an increasing 
function of the heating temperature Th and the duration of the heating process 
This reflects the fact that polymer particles will amalgamate and block the access 
of solvent to certain particles. In the extreme case we observe drug or excipient 
particles that are completely surrounded by a coat of molten polymer particles. 
This results in a certain amount of drug that is not released within a physiologically 
relevant time span. Such behavior has been observed in experimental release curves, 
see figure 7. 
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Step 3: The estimation of the distribution of path lengths, /z,(t), and the mod- 
eling of the diffusion process on these graphs is streamlined in this implementation. 
Theoretically, each molecule undergoes one dimensional Brownian motion along 
a path from its initial location until it reaches the edge of the tablet. The one- 
dimensional diffusion can be written as a limit of one-dimensional random walks, 
where the step size, Ax, and the time between steps, At, go to zero in such a 
what that, in the limit, the variances do not tend to cither zero or infinity. This 
requirement is met by requiring that 

lim v ' = 2d, 

Ax, At— >0 At 

where d is the diffusion constant, see [11, XIV. 6, p. 325]. As a coarse approxima- 
tion to the one-dimensional diffusion, therefore, we allow the molecules to perform 
random walks along the vertices of the heated contact graph G, starting at every 
vertex of type D. These random walks approximate the diffusion process that each 
drug molecule undergoes in order to exit the tablet, and therefore give us an ap- 
proximation of the distribution of exit times, fx- In practice, the physical length of 
each step in the random walk is related to the diameters of the spheres representing 
the original particles, hence the actual exit time is proportional to the sum of the 
lengths. In the simulations used to generate the distribution of exit times shown in 
figures 4 and 5, the radii of all the original particles were approximately the same, 
hence, in this case, the number of steps in the random walk was used as a surrogate 
for the time to exit the tablet. Further work must be done to investigate the scaling 
properties of these random walks, as well as the effect of including particles of very 
different sizes. However, our preliminary investigations indicate that the distribu- 
tion of the number of steps in the walks is at least qualitatively the same as the 
distribution of exit times. 

The effect of the polymer, in particular the effect of the fusing of polymer parti- 
cles in the heating process, is modeled by adjusting the probability of traversing a 
particular edge in the contact graph. Different rules can be implemented and tested. 
The simplest is to assign a conductivity Cy to the edge that connects vertices Vi and 
Vj. The conductivity is large for edges that connect vertices of type D or X, and 
small for an edge to or from a particle of type P. Let Af(i) be the set of vertices 
neighboring vertex i, Wj. Note that i S Af(i), so a molecule can remain at its initial 
position. If the walker, representing a drug molecule, is situated at vertex Vi then 
the probability of moving to the adjacent vertex Vj is given by 

p(i -► j) = dj Cij 

The walk is stopped if either an exterior vertex or a prescribed maximum number 
of steps N max is reached. Notice that the metric on G is the Euclidean distance 
between its vertices inherited from its embedding in R 3 , so that the distance a 
molecule travels on the graph corresponds to the actual distance it must travel in 
space. If a random walk does not reach one of the exterior vertices within N max 
steps, then the molecule is considered trapped inside the tablet. This will be used 
to explain partial release of drug at higher polymer concentrations. 
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3. Model 2: The continuous model. In the second model we regard time and 
space as continuous quantities and we set up a system of partial differential equa- 
tions for the contents of dissolved and undissolved excipient and drug in the tablet. 
Our starting point is the assumption that fluid quickly permeates the tablet [16] 
and that dissolution is limited by the saturation of the fluid with excipient and 
drug. This causes a delay in the initial release of drug molecules until porosity in- 
creases through diffusion at the boundary. This delay can be observed as a change 
in concavity near t = in the experimental release profiles, see figure 7. 

Let be the spatial domain of the tablet. We introduce cylindrical coordinates 
(r, 9, z) such that 

= {(r, 9,z) : < r < R, < z < H, < 6 < 2n}. 

and assume that our tablet dissolves symmetrically, i.e. concentrations do not de- 
pend on the angular variable 9 [24]. Let u± be the concentration of dissolved ex- 
cipient in the solvent. Let u 2 be the concentration of undissolved excipient in the 
solid remainder of the tablet. Likewise, we denote by v\ and v% the concentration 
of solved drug in the solvent and the content of undissolved drug in the tablet, re- 
spectively. All these quantities have dimension ML~ 3 . By k we denote the porosity 
of the tablet, k is a dimensionless value between and 1 and will increase as more 
and more excipient and drug are dissolved in the solvent. Then kui is the con- 
centration of solved excipient in the tablet. We model the diffusion classically by 
assuming that the flux of solved excipient across an interface is proportional to the 
concentration gradient times the area of the interface. The area of an interface is 
also proportional to n, i.e. the flux is given by 

-^ ua; solved excipient = ~ Du \K^Ui(r, z) ) ' 

where D u is the diffusion constant of the dissolved excipient (and D v is the diffusion 
constant of the dissolved drug) . The conservation of mass equation then yields 

d 18/ d \ d ( d \ 

(kui) = -— rD u K—Ui(r, z) + — D u n— m(r, z) + g(u 1 ,u 2 ) 



at r or \ or J oz \ Oz 

V( r , z ) ■ (D u kVui) + g(u 1 ,u 2 ), 

where g{ui, u 2 ) is the rate of concentration increase from dissolving excipient. The 
higher the porosity, the more contact there is between the solid excipient and the 
fluid in the pore space and hence the rate of dissolution increases with porosity. We 
assume this to be linear for simplicity but plan to explore more complex relation- 
ships in future work. A similar modeling assumption has been made by Lemaire et 
al. [16], see in particular equation (1) in that paper. The higher the concentration 
of dissolved excipient, the slower the rate of dissolution with no dissolution taking 
place at the saturation point C,^ ax of the fluid. Hence we assume that 

g(ux,u 2 ) = ol u k ( 1- ) u 2 (2) 

^ max 
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and we obtain the following system of evolution equations 

d_ 

Of 



^(k(u 2 , V2)u{) - V( r ^) • (D u k,{u 2 , v 2 ) Viti) = a u n{u 2 , v 2 ) (l - J M 2, 

\ nictx / 

(3) 



^-u 2 = -a u K(u 2 ,V 2 ) M - j //j 



^(k(u 2 , V 2 )V 1 ) - V( nz ) ■ {D v k{u 2 , V 2 )S/v{) = a v «(il 2 , « 2) U ~ ^ 



II If IX 



~u 2 = -a„Ac(ii 2 , v 2 ) ^1 - ^ u 2 . 

The rates of dissolution of excipient and drug are denoted by a u and a v , respectively. 
The porosity n{u 2 ,v 2 ) depends on the concentration of undissolved excipient and 
drug in the tablet. We assume that 

K(U 2 ,V 2 ) = K(U 2 + V 2 ) = (K - Kendj-fT, n + K end, 

where kq is the initial porosity and K en d is the porosity of the tablet once all the 
excipient and drug are dissolved. The initial concentrations of undissolved excipient 
and drug are denoted by u 2 and v 2 . We assume that initially no excipient and drug 
are dissolved and that the solid excipient and drug are uniformly distributed, 

u° = <u° = 0, ul{x) = u%, and v°(x) = w§, 

for positive constants u 2 and v 2 . If drug and excipient are completely undissolved 
then the initial porosity n$ will be about 2 %, that is the porosity after compaction. 
If drug and excipient are completely dissolved then the porosity will be about K en d — 
60%. Equations (3) are completed by homogeneous Dirichlet boundary conditions 
for u\ and v\ 

Ux =0, vi = 

on dfl, as we assume that any dissolved excipient or drug outside the tablet is 
immediately carried away. 

We add the first two equations of system (3), integrate over Q and apply the 
divergence theorem. We obtain, after exchanging the order of integration and dif- 
ferentiation 

— / (k(u 2 + v 2 )ui + u 2 ) dx = D u { k(u 2 + v 2 )Vui • nder =: - J u (t), 
<** Jn JdQ, 

where n denotes the outward normal and a the surface measure. This is the rate 

of change of the excipient load inside the tablet and defines the flux of excipient 

— J u (t) across the boundary d£l of the tablet. (Alternatively we could integrate the 

normal component of the flux in (1) over the boundary). A similar calculation gives 

the flux of drug 

— J v (t) = D v / k{u 2 + v 2 )Vvi ■ ndu. 
Jon 

The cumulative amount of drug released is then given by 

Rv(t) = [ J v (s)ds. 
Jo 

The solutions ui, u 2 , v\ and v 2 remain positive and 
lim (\\ Ul (t)\\ 

OO + ll u 2(£)|| + ||^2 (*) ||oo) = 0- 
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It follows then that 

lim R v (t) = \\v 2 \\i, (4) 

t — >oo 

that is, the initial drug load will be completely released. We will return to this 
point in section 5. 

4. Results. In this section we describe the implementation of our models outlined 
in sections 2 and 3. All procedures have been implemented in matlab and programs 
will be available upon request from the corresponding author. 

4.1. The Lubachevsky-Stillinger protocol. We work with a version of the 
Lubachevsky-Stillinger protocol [14, 18] that allows for classes of spheres with dif- 
ferent radii. In the first step we convert the mass fractions of drug fp>, polymer fp 
and excipient fx into particle fractions, using the known particle radii rp>, rp and 
r x and the densities Qd, Qp and q x of the pure powders. For i 6 {D, P, X} let 
m i = r i&i and set 

, - fpm D . fxm D 

n D = 1, n P = , n x = -, , 

jD-mp jDmx 

and 

rii = — — — — , 5 

n D + np + n x 

these are the particle fractions of drug, polymer and excipient in the powder mix- 
ture. The particle sizes, densities and mass fractions for different mixtures of the 
drug indomethacin (a commonly used anti-inflammatory drug), the polymer Eu- 
dragit RLPO and the excipient mannitol are given in table 1. The corresponding 
particle fractions are given in table 2. 

A difficulty in the implementation of the Lubachevsky-Stillinger protocol is the 
decision when to stop. In theory, the time between two collisions (either sphere- 
sphere or sphere-wall collisions) approaches zero and hence the increase in the pack- 
ing fraction approaches zero. We stop the execution if the times between successive 
collisions tk+i — tk, tk+2 — tk+i, ■ ■ ■ , tk+n — ifc+n-i < e, f° r a number n and a con- 
stant e > that are set by the user. The role of n is to avoid premature termination 
if accidentally a time between two collisions is very short. In our implementation we 
choose e = 10~ 6 and n = 50. This results in packing fractions of w 0.54 for packings 
of spheres with roughly equal radii (figure 2). This value is somewhat below 0.64, 
the "generally accepted" packing fraction for a random dense monodisperse packing 
[29]. However, we consider it sufficient in this preliminary work. 

At the termination of the Lubachevsky-Stillinger protocol at time t* we obtain a 
set of positions (xi(t*))fL 1 and radii (ri(t*))fL 1 . After possible compression of this 
sphere packing, we define the graph G by joining vertices Xi and Xj if 

\Xi -Xj\ < A(r. i + r J ) (6) 

where A > 1 is a constant. Graphs of this kind of graph are also known as proximity 
graphs on random point sets [19,20]. The heating process removes edges with a 
probability p and results in the heated contact graph G. Figure 3 shows a contact 
graph of a packing before and after heating. Choosing A slightly greater than 1 
takes into account that the molecules can travel through interstitial spaces between 
particles that are close, but not touching. The spheres that have at least one surface 
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point close to a boundary of the domain are given the "exterior" label. For example, 
if 

x x 3 < firj or > 1 - fir 3 (7) 
for a constant /i > 1, then that sphere is a possible point of exit in the x-direction. 



D:P:X = 90:100:545 D:P:X = 90:100:545 




Figure 3. The proximity graphs on the sphere centers of a ran- 
dom dense sphere packing. The constants in equations (6) and (7) 
are chosen A = /i = 1.03. The left panel shows the contact graph G 
before the heating process. The edges are removed with a proba- 
bility p = 0.3 during the heating process, this results in the heated 
contact graph G on the right. 
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drug 


polymer 


excipicnt 


diameter (/xm) 


100 


100 


106 


density (g/cm 3 ) 


1.37 


1.23 


1.52 


diffusion coefficient (cm 2 / s) 


7- 10~ 6 




7- 10~ 6 



Table 1. Composition of tablets, particle sizes and true densi- 
ties for drug (indomethacin) , polymer, and excipient (mannitol or 
lactose). The fractions are the mass fractions in % in a tablet of 
total weight 500 mg. The heights hi and hi are those of a tablet 
of 13 mm diameter at compression pressures of 74 M Pa (hi) and 
221 MPa(hi), respectively. The diffusion coefficients of drug and 
excipient are those in water. 
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n D 


Hp 


n x 


12.2 


13.6 


74.2 


11.7 


26.1 


62.2 


11.3 


37.6 


51.1 


10.8 


48.2 


41.0 


10.4 


58.0 


31.6 



Table 2. The particle fractions of drug, polymer and excipient (in 
%) corresponding to the mass fraction in table 1 are computed as 
in equation (5). 



4.2. Release curves obtained from random walks on the contact graph. 

On the heated contact graph G we perform random walks that start from each drug 
particle and end when an exterior vertex is reached, or the maximum number of 
steps, N max , is achieved. As explained above, the probability of traversing edge 
is determined by the type of the terminal vertex, 

f c_ if j is of type P, , g . 
13 1 c+ otherwise 

we use here c + = 100 and c_ = 1. Notice that the current position is also a possible 
terminal vertex, that is, the random walker can remain in a fixed place for (perhaps 
long) periods of time. For every path we record the number of steps in the path 
and the length of the edges traversed. If a particle remains at a vertex, we still 
add "one" to the number of steps, and we add the diameter of the current vertex' 
particle type to the path length. Thus, even though a molecule may be close to 
where it started after one time step, the time step is still counted as "time spent 
before exiting" . 

We present simulated release profiles in Figures 4 and 5. In these figures we vary 
the polymer fraction and keep all other parameters constant. We observe that as 
the polymer fraction increases, the drug release is slowed down. Further, we test 
different edge removal probabilities (p = 0.3 and p = 0.5) that can be interpreted 
as different curing temperatures. An increase in p decreases the fraction of drug 
released. 

4.3. Release curves predicted by the continuous model. As a first simpli- 
fication we disregard the height of the tablet and collapse it to a disk, so that 
functions are now only dependent on the variable r £ [0, 1]. We split the numer- 
ical solution procedure into a diffusion step for U\ and v\ and a reaction step for 
all four concentrations. It is assumed that during the diffusion step the porosity 
k(u 2 + V2) does not change. We discretize the radial Laplace operator on the uni- 
form grid rk — (k— l)Ar, k — 1, . . . , N using standard finite differences and we use 
the Crank-Nicolson scheme at every diffusion step [23, Section 17.3]. The reaction 
step is calculated using the Euler forward method. Release curves obtained from 
numerical solution of equation (3) are shown in figure 6. A larger final porosity 
results in a faster release of the drug. Observe that all release curves in figure 6 
predict a complete release of the drug, leading to the conjecture (4). 

4.4. Comparison with experimental release curves. Experimental tablets 
were formulated from mixtures of indomethacin, Eudragit RLPO and mannitol 
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Figure 4. Five release profiles obtained from step numbers of ran- 
dom walks as described in section 4.2, as the polymer mass fraction 
varies. The mass fraction of drug is 10% throughout and their num- 
ber in each packing is 180. For the actual particle fractions of drug, 
polymer and excipient see the left part of table 2. We have used 
A = [i = 1.03 in equations (6) and (7) and a probability for edge 
removal p = 0.1. The edge propensities are given in equation (8). 
The maximum number of steps is N max — 10 3 . 
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Figure 5. As in figure 4, but now with a probability for edge 
removal p = 0.3. 



and their release profiles were determined, see [8] for a detailed description. Briefly, 
powder mixtures were compressed at different pressures and tablets were heated at 
selected temperatures. The tablets were then placed in a phosphate buffer medium 
and samples were collected at different time points over a period of 8 h (see figure 
7) . We observe a good qualitative agreement with our simulated release profiles in 
figures 4 and 5. 

5. Discussion and Conclusion. Several studies have investigated percolation on 
regular lattices and graphs obtained from random sphere packings. Villalobos and 
collaborators [30,31] used Monte Carlo simulations on cubic lattices to predict drug 
release profiles from binary drug/excipient mixtures. In these works the excipient 
played the role of the inert matrix, which in our case is the polymer. It was re- 
ported that if the particle fraction of the inert matrix (the inaccessible sites in the 
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Figure 6. Release profiles predicted by the continuous model (3) 
as the final porosity K en d varies. The dimensionless parameters 
used in this example are k — 0.02, D u = 0.3, D v = 0.5 a u = a v — 
1.5 and C^ ax = C^ ax = u 2 = v 2 = 1. 



lattice) is below a threshold of w 69%, then the release of the drug is slowed down 
with increasing matrix fraction, but will still be complete. Only above the critical 
matrix concentration an entrapment of drug will be observed. The value of 69% is 
complementary to the estimated site percolation threshold of 0.3116 for the cubic 
lattice [26]. A similar value was reported by Powell [21,22] for percolation on the 
contact graph of a monodisperse random sphere packing. In this work we include 
the heating of the tablet, which melts the polymer. The blocked polymer particles 
and the removal of certain edges provide a case of what is known as site-bond perco- 
lation [7] . In their paper [7] , Chang and Odagaki studied site and bond percolation 
processes where sites and bonds are removed independently from each other and 
eventually open sites are connected if they share an open bond. They found that 
for the cubic lattice, the percolation threshold as a function of the probabilities 
that sites respectively bonds are open, is well described by a hyperbola. While 
this was not the primary goal of our study, we ended up addressing the problem 
of determining percolation thresholds for contact graphs of dense random sphere 
packings. 

While we have seen a good qualitative agreement between our simulated release 
profiles and experimental data, further research is needed to provide more reliable 
and quantitative predictions. In our simulations we find many particles with very 
short escape paths. However, their number is going to diminish as we consider 
sphere packings with more spheres. Intuitively, the fraction of spheres "close to the 
exterior" is decreasing as the number of spheres in the packing increases, a problem 
that was discussed already in [21]. Larger simulations, perhaps using parallel com- 
puting methods, are required. Further work will also include careful calibration of 
the parameters A, /i and p to experimental release curves. 

Our continuous model captures nicely the change of the release curves from con- 
vex to concave. However, as figure 6 and equation (4) suggest, the release will 
always be complete. No percolation behavior is exhibited by the system (3). In 
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Figure 7. Release of indomethacin (mass fraction 10%) from Eu- 
dragit RLPO matrix tablets containing mannitol (90 — 125 /im par- 
ticle diameter), a plastic excipient, using the USP basket apparatus 
at 100 rpm in 900 ml phosphate buffer pH = 7.2 (0.2 M) at 37° C. 
The top figure shows a tablet that was compressed at 74 MP a and 
cured 24 h at 40° C, the bottom figure shows a tablet that was 
cured 24 h at 70° C. Data are representative for three repetitions 
of the experiment (from [8]). 

the future we need to take into account the permeability of the porous matrix, a 
notion that arises in hydrology [2]. In the derivation of the dissolution kinetics 
(2), we have assumed that the dissolution of drug and excipient is helped by an 
increase in porosity. This is contrary to the commonly made assumption of a reced- 
ing boundary in the pharmaceutical literature [12]. There, the drug dissolves in a 
way that decreases the area of contact between drug and solvent. It poses an inter- 
esting inverse problem to distinguish these two concepts, given the experimentally 
determined release profiles. 

Another interesting suggestion of Villalobos et al. is that the release profiles have 
a Weibull cumulative distribution function, [30, Equation 2]. 



^ = 1 - cxp(-at b ) 

OO 
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where the parameters are influenced by the size of the matrix tablet. More precisely, 
the ratio of the number of exit sites versus the total number of sites, N exit /N should 
influence the release kinetics, with faster release for smaller values of N ex u/N. If 
we assume that N exit oc N$ then larger tablets result in slower release. While this 
makes sense intuitively, we consider this to be an interesting question to explore 
more rigorously, and we plan to address it in future work. 
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